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Abstract 

This paper outlines the development of a mathe- 
matical model that is expected to be useful for ro- 
torcraft flying qualities research. A computer model 
is presented that can be applied to a range of dif- 
ferent rotorcraft configurations. The algorithm com- 
putes vehicle trim and a linear state-space model of 
the aircraft.. The trim algorithm uses non linear op- 
timization theory to solve the non linear algebraic 
trim equations. The linear aircraft equations consist 
of an airframe model and a flight control system dy- 
namic model. The airframe model includes coupled 
rotor and fuselage rigid body dynamics and aerody- 
namics. The aerodynamic model for the rotors uti- 
lizes blade element theory and a three state dynamic 
inflow model. Aerodynamics of the fuselage and fuse- 
lage empennages are included. The linear state-space 
description for the flight control system is developed 
using standard block diagram data. 

Introduction 

In the past, rotorcraft flight control system pre- 
liminary design used mathematical models which as- 
sumed the fuselage to possess six degrees of freedom. 
The rotor dynamics were assumed to be substantially 
faster than the fuselage dynamics and were subse- 
quently approximated as quasi-static. The process 
of fine tuning the flight control system was accom- 
plished through an extensive flight test program com- 
prised of a matrix of control system parameter varia- 
tions. While fine tuning of the flight control system is 
still accomplished through flight testing the vehicle, 
significant improvements in the optimization process 
have been realized when high order dynamic rotor- 
craft models are utilized during the preliminary flight 
control system design stage. 

•Presented at Piloting Vertical Flight Aircraft: A Confer- 
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Rotorcraft are now being designed with sophisti- 
cated electronic flight control systems. These com- 
plex control systems are utilized not only to satisfy 
standard flying qualities specifications but also to 
meet aerodynamic performance, vibration, and struc- 
tural loads criteria. The design of modern rotorcraft 
flight control systems now stretches across many dif- 
ferent individual disciplines and is indeed interdisci- 
plinary. The general trend toward increased reliance 
on the flight control system for improving overall sys- 
tem performance has lead designers to consider higher 
bandwidth systems which rely on high levels of sensor 
feedback to yield desired aircraft stability. The main 
drawback of this approach is that increased levels of 
feedback, which in general improve the low frequency 
fuselage dynamic behavior, can destabilize higher fre- 
quency rotor blade motion. In order to make mean- 
ingful estimates of the impact of a particular flight 
control configuration on system requirements it has 
been found that a mathematical model which in- 
cludes fuselage and rotor rigid body dynamics and 
rotor dynamic inflow is necessary [1]. 

The business of rotorcraft modeling for flight con- 
trol system design and analysis support has been an 
active research area for many years. Deriving the 
equations of motion of a fully coupled fuselage and 
rotor system for a reasonably general configuration 
quickly becomes unwieldy due to complicated ge- 
ometry including many matrix transformations and 
intricate logic branching. These complexities have 
lead engineers to develop digital computer programs 
which more or less relegate model computation to the 
computer and free the engineer to focus on analysis 
results. 

Talbot, Tinling, Decker, and Chen [2] formulated 
a helicopter flying qualities model that includes fuse- 
lage dynamics and a three degree of freedom tip-path- 
plane representation for the main rotor flapping dy- 
namics. Some simplifications are made in the anal- 
ysis in order to formulate compact, analytical force 
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and moment expressions for the rotor forces and mo- 
ments. Gibbons and Done [3] derived a numerical 
method to automatically generate rotorcraft equa- 
tions of motion. The method uses Lagrange’s equa- 
tions and relies on expressing inertial position vectors 
of the rotor blades as a matrix multiplied by the po- 
sition vector in blade coordinates plus a term that is 
a function of the modal coordinates, time, and span- 
wise position. The required differentiations of the 
position vector to form the equations of motion are 
performed numerically. Miller and White [1] used 
concepts from Lytwyn [4] and Gibbons and Done 
[3] to automate generation of the equations of mo- 
tion for rotorcraft handling qualities analysis. Miller 
and White [1] expressed all transformation matrices 
in complex variable form and were able to develop a 
compact algorithm to analytically obtain long strings 
of orthogonal transformation matrices along with all 
necessary derivatives to form nonlinear and linearized 
dynamic equations. Lagrange’s equations were used 
in the formulation. Zhao and Curtiss [5] derived a 
set of linearized equations by analytic linearization of 
a nonlinear model formulated using Lagrange’s equa- 
tions. The symbolic manipulation computer program 
MACSYMA was used in forming the equations. Sub- 
sequent work by McKillip and Curtiss [6] has im- 
proved and extended the work by Zhao and Curtiss 

[ 5 ]- 

The work discussed in this paper derives a rotor- 
craft flying qualities model which has been imple- 
mented into a FORTRAN computer program. A 
fairly generic rotorcraft configuration, consisting of 
a rigid fuselage, two rotors, and an arbitrary number 
of fuselage fixed external surfaces has been assumed, 
as shown in Figure 1. It is important to note that the 
type of analysis carried out in this work can accom- 
modate any arbitrary number of rotors in the con- 
figuration. The number of rotors has been chosen to 
be two since the majority of rotorcraft fall under this 
category. The fuselage possesses six degrees of free- 
dom and the rotor blades have flap, lag, and pitch 
degrees of freedom. The rotor aerodynamic models 
are based on blade element theory and include three 
degree of freedom dynamic inflow. The equations 
of motion are formulated using Kane’s equations [7]. 
More importantly, derivatives of transformation ma- 
trices are formed using angular velocity expressions 
as opposed to numerical or direct differentiation. The 
rotor dynamic inflow equations are based on the Pitt 
and Peters model [8] and include hub motion pertur- 
bations. The residual of the equations of motion and 
the residual gradient expressions are derived analyt- 
ically and trim is calculated using the residual and 
residual gradients in concert with a modified New- 


ton’s method. The rotor trim variables are the ro- 
tor multiblade coordinates. A linear constant coeffi- 
cient model of the composite airframe is formulated 
using a multiblade coordinate transformation with a 
subsequent constant coefficient approximation. The 
linear constant coefficient airframe model is coupled 
to the linear control system dynamic model to form 
the overall linear model. Linear analysis tools such 
as eigen values, eigen vectors, transfer functions, fre- 
quency response, and linear simulation are directly 
contained within the computer program. 


Airframe Dynamic Model 

As pictured in Figure 1, the airframe dynamic 
model consists of a rigid fuselage with the standard 
six degrees of freedom and two fully articulated ro- 
tor systems, each with dynamic inflow. The fuse- 
lage aerodynamic force and moment components are 
obtained in the wind axis from a two dimensional 
data table as functions of fuselage angle of attack 
and sideslip. The aerodynamic forces exerted on the 
external surfaces are obtained using standard lifting 
line theory. The rotor geometry details are shown in 
Figure 2. Provisions are made in the model to accom- 
modate any of the six possible sequences of flap, lag, 
and pitch hinges for the rotor blades. Each hinge is 
accompanied by a linear torsional spring and damper. 
Each blade also has a non linear translational damper 
which is attached to the rotor blade from the rotor 
hub. Hingeless rotor systems can be approximately 
modeled using a virtual hinge representation. The 
aerodynamic forces exerted on the rotor blades are 
calculated using blade element theory. The blades 
on a rotor have identical yet arbitrary geometric and 
inertial properties. 

The airframe nonlinear dynamic model is obtained 
using the flat and non-rotating earth assumption. 
Kane’s Equations are then written for each degree of 
freedom by taking into account the contributions of 
the generalized inertia forces, the generalized gravity 
forces, the generalized aerodynamic forces, and the 
generalized spring-damper forces. 

/r(0 = fl r (t) + f Gr (t) + f Ar (t) + f SDr (t) } 

r = 1, . . . , n RB (1) 

In equation 1, t denotes time and n RB is the number 
of generalized speeds. The origin of each term on the 
right hand side of Equation 1 is discussed below. 

The following nomenclature is introduced for deriv- 
ing the generalized inertia forces. Let n R1 and n R2 
denote the number of blades on rotor 1 and rotor 
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2, respectively. Let m/p and /f, m/h.i and Irij (i = 
1, . . . , n R i), and m R2tj and I R2J (j - 1, — , n^), re- 
spectively, denote the masses and inertia matrices for 
the fuselage, rotor 1 blades, and rotor 2 blades. Let 
wf, URi ti (i = 1, . . . , Tim), and u R2 j (j = 1, — , n/? 2 ) 
represent the individual body axis components of the 
angular velocities of the fuselage, rotor 1 blades, and 
rotor 2 blades, respectively. Let and a/p., 
and (* = l,...,n/*i), and v R 2 j* and a R 2 j* 

(j — 1, ^,ti R2 ) represent the inertial axis compo- 
nents of the c.g. (center of gravity) velocities and ac- 
celerations of the fuselage, rotor 1 blades, and rotor 
2 blades, respectively. Then the generalized inertia 
forces acting on the configuration can be written as, 
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{I R2 R2 ,i + S(Wfl2,i)Jfl2 l * w rt2,t} , 


r = 1, . . . , n RB (2) 


where an overdot denotes differentiation with respect 
to time and 5(*) is the standard cross product skew- 
symmetric matrix operator (Appendix), u is the vec- 
tor of generalized speeds. Letting g be the accelera- 
tion due to gravity, the generalized gravity forces can 
be written as, 


fa r (t) 


-m F g 


9(vf * )3 
du r 


kri 

~E mfl M g 
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djvR i,i»)3 
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r = 1, . . .,tirb (3) 

The generalized aerodynamic forces are discussed 
next. Let vp, and F f and M F be respectively 
the body axis components of the velocity and the 


aerodynamic force and moment acting on the fuse- 
lage aerodynamic center. Let b R2 , and bs, 

( i = 1 ,...,715) denote the number of elements or 
sections on any rotor 1 blade, any rotor 2 blade, 
and the ith external surface. Let v R \ t ij and Fju.tj 
(i = - l,... 1 6iu), and v R 2 l k,i and 

FR 2 t h t i (i = 1 ) • • • j n R 2 , J = l>-"i(A 2) be the indi- 
vidual body axis components of the section velocity 
and the aerodynamic force acting on rotor 1 blades 
and rotor 2 blades, respectively. Let vs t j and Fs t j 
(i = 1 , . . . , 7 %s,j = 1 , . . . , bs t ) be the respective body 
axis components of the section velocity and the aero- 
dynamic force acting on the external surfaces. The 
generalized aerodynamic forces acting on the config- 
uration can then be written as, 



r = 1, . . . , tirb (4) 

The generalized spring-damper forces are discussed 
next. Figure 2 shows the typical spring-damper at- 
tachment geometry for a typical blade. Let v Rl • • 
and F Ri,ij (*' = J = 1,2), and vg 2tl 

and Ff( 2kl ( k = 1,... , ti R 2 j f = 1,2) be the indi- 
vidual rotor hub axis components of the velocities 
and the forces acting on the translational damper 
attachment points for rotor 1 blades and rotor 2 
blades, respectively. Let u R nj and M Rli j (i = 
= 1 ,..., 4 ), and w% 2kt and M% 2il 
(k = = l,-.., 4) denote the individual 

body axis components of the angular velocities and 
the torsional spring-damper moments acting on the 
hub, link 1, link 2, and the blade for rotor 1 blades 
and rotor 2 blades, respectively. Then the generalized 
spring-damper forces can be expressed as, 


fSDr(t) = “EE 


nR1 2 f fhp 

^ / av Rl,iJ 

du r 
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r — 1 ,---,n R B (5) 

The partial derivatives and in Equations 
2 through 5 are known as partial velocities and par- 
tial angular velocities, respectively. The generalized 
coordinate vector q and the generalized speed vector 
u are defined as follows: 

(?R2,l) T , (9ft2,2) T , • • •,(?R2,n R2 ) T } T (6) 

“ = {( U f) T .( U «l,l) T >( U fll,2) T ,...,(«fll. nRI ) T , 

(wR2,l) T , («R2,2) T . • • • , (uR2, nR2 ) T } T (7) 

The subscripts F } Rl , and R 2 refer to fuselage, rotor 
1, and rotor 2 variables, respectively. Further, 
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fuselage axes, T H \,i ( i = is the matrix 

of transformations from rotating hub axes to shaft 
axes, Tfti t (i = 1 , . . . , tiri) is the matrix of transfor- 
mations from link 1 axes to rotating hub axes, 

(i = 1 , . . . , n/ji) is the matrix of transformations from 

link 2 axes to link 1 axes, and t£ 3) • (i ~ 1, n m ) is 

the matrix of transformations from blade axes to link 
2 axes. The transformation matrices are expressed as 
follows: 


T f 

= [E 1 (<t>)E 2 ( 6 )E 3 (rP)f 

(14) 

T S 1 
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In Equation 16, rp Rl>i = Q m t + ^(i - 1), where 
Qri is the rotor 1 hub rotational speed. It is as- 
sumed that the shaft is inclined with respect to the 
fuselage by first a rotation with the angle r$i )a and 
then a rotation with the angle rsi,*. Depending on 
the sequence of rotation, T S \, a and T SX)b are one each 
among E\ and E 2. E Xl £2? and £3 are single axis 
transformation matrices about x, y , and z axes, re- 
spectively (Appendix). Clearly, T^\ T&\ anc j ji$) 
are one each among E\ y E 2, and £3, depending on 
the rotor blade hinge sequence. 

The body axis components of the angular velocity 
of the fuselage, can be written as, 


The quantities a< 1 ), a (2) , and a & ar e one of lag, flap, 
and pitch angles, depending on the rotor blade hinge 
sequence. 

A brief description of the analysis involved in calcu- 
lating the terms on the right hand sides of Equations 
2 through 5 is given in the following. For simplicity, 
the analysis for the rotor terms will be restricted to 
rotor 1; the analysis for rotor 2 terms is analogous. 

Generalized Inertia Forces 

The six terms comprising the generalized inertia 
forces, Equation 2, are discussed here. The orienta- 
tion of the fuselage with respect to inertial axes and 
the orientation of the rotors with respect to the fuse- 
lage can be described using transformation matrices. 
Each transformation matrix is composed of one, two, 
or three single axis transformation matrices. Refer- 
ring to Figure 2, let Tp be the matrix of transfor- 
mations from the fuselage axes to inertial axes. The 
following five matrices are defined for rotor 1. Ts\ 
is the matrix of transformations from shaft axes to 


UF = {p, q,r} 


(2U) 


Using the transformation matrices defined above, the 
body axis components of the angular velocity of the 
ith rotor 1 blade can be written as, 




LJF + 


[osy : 


°Rl a Rl,i + 


[jt( 3) 1 T r(2) • (2) 

°Rl a Rl t i + 


+ 


& (3 >d (3) / 0n 

°Rl a Rl,i (21) 

The unit vectors 6^, and b ^ have been intro- 
duced to allow a general rotor blade hinge sequence. 
For example, if rotor 1 blades undergo a lag, flap, 
and pitch rotation sequence, then 6^ = {0,0, 1} T , 
^ri — and b ^ = {1,0, 0} T . Equation 21 
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has been obtained using the concept of simple angular 
velocities [7]. 

The body axis components of the fuselage c.g. ve- 
locity are given as, 


v b = {u,v,w} t ( 22 ) 

The inertial axis components of the fuselage and 
blade c.g. velocities can be written as, 

Vp* — TpVB ( 23 ) 

VRl } i * = vp * + 

TpS(u F )rpi + 

TFTs\S(wsi)fHl + 
T F TsiT ffU S(u Hiii )r ( £ -f 
TFTsiTHuTgliS^l)^ + 

T F TsiT Hl , i T^ 4 T<£ 4 S(w^ i yS + 

TfTsiThu^tS^S^M, 


(24) 


In Equation 24, u>si represents the body axis compo- 
nents of the angular velocity of rotor 1 shaft. 

and ■ are the individual body axis com- 
ponents of the angular velocities of the rotating hub, 
link 1, and link 2, respectively. The expressions for 
these angular velocities are given as follows: 
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In Equation 24, the vectors rp\ t tr 1 , r^, 
and are defined as follows. fp\ is the position vec- 
tor from fuselage c.g. to a point on shaft 1, expressed 
in fuselage axes, tr\ is the position vector from the 
point on shaft 1 to the center of hub 1, expressed 
in shaft 1 axes. For any rotor 1 blade, is the 
position vector from the center of hub 1 to the first 
hinge, expressed in rotating hub 1 axes; is the po- 
sition vector from the first hinge to the second hinge, 
expressed in link 1 axes; is the position vector 
from the second hinge to the third hinge, expressed 
in link 2 axes; and is the position vector from the 
third hinge to the blade c.g., expressed in blade axes. 
Equations 20 through 24 are used to compute the par- 
tial velocities and partial angular velocities needed in 
Equation 2. 


The angular acceleration vectors ujp and u; dis- 
appearing in Equation 2 are obtained by a time- 
differentiation of the right hand sides of Equations 20 
and 21, respectively. Similarly, the translational ac- 
celeration vectors ap • and appearing in Equa- 

tion 2 are obtained by a time-differentiation of the 
right hand sides of Equations 23 and 24, respectively. 
While the equations for the rotor blade acceleration 
vectors are lengthy and omitted here, it is noticed 
from an inspection of Equations 20 through 24 that 
obtaining these equations is straight forward once the 
expressions for the time-derivatives of the transfor- 
mation matrices has been obtained. The Appendix 
gives the derivation of a formula for calculating the 
time-derivative of a matrix in terms of a matrix prod- 
uct. Using this formula, the following are obtained: 
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= TfS(uf) 

(29) 

Tsi 

= 0 

(30) 

Tm,i 

= T H i,iS({ 0,0,-Q fll } T ) 

(31) 
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— \°Ri a Ri t i J 
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ft 2) 
i m,i 

— T ^ S ( ^ 

\°Rl a Rl t i ) 
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r( 3 ) 

1 Rl,i 

— T ^ S ( ^ 

“ 1 Rl,i^ { 0 Rl a Rl t i ) 

(34) 


Generalized Gravity Forces 

The partial velocities 1 an< ^ ob- 

tained in the computation of generalized inertia forces 
are used to compute the generalized gravity forces 
given by Equation 3. 

Generalized Aerodynamic Forces Due to Fuse- 
lage 

The first term in Equation 4 represents the gener- 
alized aerodynamic forces due to the fuselage. The 
quantities comprising this term are obtained as fol- 
lows. The body axis components of the velocity of the 
fuselage aerodynamic center (a.c.) can be written as, 

v f = v b + S(u>p)r A c (35) 

where f A c is the position vector from the fuselage c.g. 
to the fuselage a.c., expressed in body axes. Equation 
35 is used to compute the partial velocity . 

For a rotorcraft, the wind- axis components of the 
aerodynamic force and moment acting at the fuselage 
a.c. are usually given as a function of fuselage angles 
of attack and sideslip: 


L = 

Li(a) + L 2 (0) 

(36) 

D = 

£>i(«) + W) 

(37) 

M = 

M 1 (a) + M 2 (/3) 

(38) 
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y = Yi(q)+y 2 (0) 

(39) 

> = h(a)+h(0) 

(40) 

N = N!(a) + N 2 (0) 

(41) 


These forces and moments are scaled with respect to 
the local dynamic pressure and can be in the form 
of a two dimensional data table or fitted analytical 
expressions to wind-tunnel data. The force and mo- 
ment components in the body axes are given as, 


T S( US i)r H i + 

[ 4 i!, 4 ?,] T 5 ( 4 V,Jr^ + 

S(uRi,i)f RlJ (4g) 


M f = qE 2 (a)E 3 (-/?) {/, M, N f 

The fuselage velocities, for purposes of calculating the 
aerodynamic variables a, 0, and q, include the effect 
of rotor 1 down wash: 


F f = qE 2 (a)E 3 (-0) {-D,Y,-L} t (42) The ° nly ne T 9 uant >ty introduced in the preceding 

Equation is r R1 j (j — l,...,6 fil ), which is the po- 
(43) sition vector from the root of the blade to the ;'th 
aerodynamic element, expressed in blade body axes. 
f^ at ' 0n 4® * s use< i 1° obtain the partial velocity 

0u r 1 ' 

Figure 3 shows a typical jth element on the ith 
blade, and the lift and drag forces acting on it. 
(y«i,i, are the body axes of the blade. 9 Rl j 

is the blade twist angle at the ;'th section, , 

and u Ri,ij are components of the relative air ve- 
locity parallel and perpendicular to the zero lift line. 
The variables a Miijl and D Rl ,i j have ob- 

vious meanings. The velocity of the (i,j)th element 
with respect to air is given by the following equation, 


v — v F — in«i,o/fii(x/?i) 


(44) 


where W } ji q is the rotor 1 collective inflow, and 
is the rotor 1 wake skew angle. In absence of more 
sophisticated data, f R1 assumes the value (0, 0, 1} T 
or {0, 0, 0} , depending on whether the a.c. is within 
or outside of rotor 1 wake, y/jj is given as, 


X/u = tan -1 


Pfii 


-A 


R 1 


(45) 


where fi m and X Rl are, respectively, the rotor 1 ad- 
vance ratio and rotor 1 inflow ratio. These quantities 
can be determined by computing the relative air ve- 
locity components at the rotor hub. Using Equation 
44, the aerodynamic variables a, 0, and q can be 
readily computed: 


a — tan 


“(if)' 


0 = sin 1 ( 7 ^~ 
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q= 2 p 


(46) 


(47) 


= vruj + [T ff i44V,,^4 3 M] T ^i.,. / 

rpi . . J . (49) 

ine term v R1 arises due to rotor inflow and can be 
approximately evaluated as, 

( t, Hl,ij)l = 0 

( v Rl,i,j)2 = 0 

(®jti,ij)3 = -%i,o-[Am+(fflij)i/fl fil ]. 

(u'fli.ijsin V’m.i + wri,i c cosiI> Ri ,) 

where (50) 


_ ( f Rl)l + (r ( R1 


hRl = 


ai)i + (yff)i 
Rr 1 


(51) 


Generalized Aerodynamic Forces Due to Ro- 
tors 

The second term in Equation 4 represents the gen- 
eralized aerodynamic forces due to rotor 1. Blade ele- 
ment analysis is used to calculate this term. As men- 
tioned earlier, the rotor blades are allowed to have 
any arbitrary variation of twist, chord length, and 
airfoil characteristics along the span. The body axis 
components of the blade element velocity are given 
as, 


wm.ij and are the sin and cos components of 

the rotor inflow and R R1 is the rotor radius. The 
radial, tangential, and perpendicular components of 
the air velocity at the airfoil can be computed as, 

{“fll.ij! -u fll,iji“/ti,, j} r = El (0Rl,j)vRl,iJ (52) 

Using the air velocity components, the section an- 
gle of attack and Mach Number can be calculated as 
follows: 



[TrTsiTwTigrfSsTjg, 

\T 

a Rl,ij 
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l Tfj i t i Tft ^ J S(u F )r pi + 
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where c is the speed of sound at the altitude where 
the aircraft is operating. Airfoil lift and drag coeffi- 
cients are usually specified as a function of the angle 
of attack and Mach Number. Thus, 

= cJ m y i tJ ( a Rl } i,jy^Rl,ij) (55) 

c m,«j = CRi,ij( a Ri,ij>MRi,ij) (56) 

The above data can be either in the form of a two 
dimensional data table or in the form of fitted an- 
alytical expressions to experimental data. However, 
in the absence of any data, simple analytical lift and 
drag models can be used. Based on reference [9], 
equations were generated for two simple models, one 
that ignored stall and compressibility effects and an- 
other that included the same. The section lift and 
drag forces are computed next: 

Dm,i,j = tfRi.ijCfli ,ij c Ri (58) 

where 

= jpKufcuj)* + («fii, t ,j) 2 ] (59) 

and CRij and tjjh j are, respectively, the chord length 
and lift efficiency factors at the jth section. The body 
axis components of the section aerodynamic force are 
given by: 

(60) 

Generalized Aerodynamic Forces Due to Sur- 
faces 

The generalized aerodynamic forces due to the ex- 
ternal surfaces are derived in much the same way as 
those due to the rotors. One difference, however, is 
that the radial drag force due to radial flow is consid- 
ered here. It is assumed that every surface has a fixed 
(invariant with time) orientation with respect to the 
fuselage. The orientation can be specified uniquely 
in terms of rotation angles about three mutually per- 
pendicular axes. In order to have consistency in de- 
scribing surfaces with different orientations, the body 
axes for any surface are defined as follows. The x axis 
coincides with the zero lift line of the root section 
and is directed from the trailing edge of the surface 
to the leading edge of the surface (see Figure 4). The 
y axis is perpendicular to the x axis, passes through 
the aerodynamic center of the root section and is di- 
rected outboard, z axis completes the right handed 
set. 

The objective is to evaluate the last term in Equa- 
tion 4. For illustration purposes, the mathematical 


analysis involved is outlined for surface 1. Let the 
surface be oriented with respect to the fuselage by 
three successive rotations of angles js x , 5s x and es x 
about mutually perpendicular axes. The sequence of 
rotations can be any of the possible six sequences. Let 
the matrices associated with the above transforma- 
tions be and Tc?\ respectively. The ma- 

trix T E5l = transforms 

components of a vector from fuselage axes to surface 
I axes. Let rs l be the position vector from the fuse- 
lage c.g. to the surface I reference point, expressed 
in fuselage axes. Let rs lt j be the position vector from 
the surface reference point to the aerodynamic cen- 
ter of the j th section, expressed in surface axes. Then 
the velocity of the jth section can be expressed as, 

v s lt j = TeSx [vb + S(up)r Sl + 5(u;f)[Te5 1 ] T 7 ; 5 1 j] 

( 61 ) 

This expression is used to obtain the partial velocity 
dV Q „ 1 ' 1 , which is needed for evaluating the generalized 
aerodynamic force contribution from surface 1. 

For purposes of computing the aerodynamic force, 
the resultant section velocity with respect to air in- 
cludes the effect of rotor 1 downwash: 

VS,J - v S lt j - TES^RlfifsAXRl) ( 62 ) 

Similar to the case of the fuselage, in a simple anal- 
ysis, f$ x can be taken to be equal to {0,0, 1} T or 
{0, 0, 0} T , depending on whether the surface is within 
or outside of rotor 1 wake. 

The effect of radial flow on a surface section is in- 
cluded in the same way as described in reference [10], 
where profile power is computed due to radial flow 
at blade sections. Let the free stream velocity at the 
j th section be yawed, as shown in Figure 5. An esti- 
mate of the normal and radial drag forces is desired, 
preferably in terms of the two dimensional sectional 
aerodynamic coefficients. It is assumed that the total 
viscous drag on the yawed section acts in the same di- 
rection as the free stream velocity. It is also assumed 
that the yawed section drag coefficient is given by the 
two dimensional unyawed airfoil characteristics. The 
normal section lift coefficient is assumed not to be 
influenced by yawed flow. The angle of attack and 
Mach Number for the unyawed and yawed sections 
are, 



= tan 1 ^ 

.fe,j)i7 

(63) 

Ms u j 

\A«Sij 

c 

(64) 


= tan -1 ^ 

( S S,j) 3 \ 

V(»S.j)l + (®5,j)5/ 

I (65) 
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m 3iJ = (66) 

The section lift and drag coefficients are given by: 

c Sij = c 5ij( a 5ijiA<Sij) (67) 


is the position vector from the blade root to attach- 
ment point 2, expressed in blade axes. The preceding 
equation is used to determine the velocity of attach- 
ment point 2 relative to that of attachment point 1, 
expressed in hub axes: 


c 5i,j - c s l ,j( < xs u j,Ms i ,j) (68) 

As mentioned for the case of rotor aerodynamics, in 

the absence of lift and drag coefficient data, simple 
analytical models for the coefficients can be used. 
The section lift and drag forces are given as, 

,i = fo,jC5 1J -Cs 1 j(Af5 l j)2^5 ll j (69) 

J ~ Qs l ,jC d s u jC Sl ,j(Ar Sl ,j )2 (70) 

where 

9Si j = j/>[(®s,j)* + (vs,j)i] (71) 

9S*J = 5 /»l(®s,j)? + («5 t j)| + («s t j)|] (72) 


VRl,i 




+ 


T Rl,Tk%T ( R %S(u R1 /)t Rl 


(75) 


The component of this relative velocity along the 
damper arm can be written as. 


VRi,i — (1/ j dfu t i \)(dfU i) T v Rl i (76) 

The damper force Frij is assumed to be specified as 
a function of the above speed. Hence, 


and rjSij is the section lift efficiency factor. Finally, 
the body axis components of the section aerodynamic 
force are given as, 


R Sui = (l Si ,i/yj (®s, ~j )? + (tS Sl ~j )l \ ■ 

{( c s t ,i) 3, 0, — (v Sl> >)i} 7 ' + 

{°Si j /yj (®s, j )i + (vs, ,j)l + (t/ s , J )§) • 
{~(®5ij)li “(®Jlj)2, -(®5,,i)3} T (73) 


Generalized Damping Forces Due to Transla- 
tional Dampers 

As shown in Figure 2, one end of the blade trans- 
lational damper is attached to the rotating hub while 
the other end is attached to the blade itself. The 
damper force is assumed to be given as a function of 
the relative speed between it’s two ends. The analysis 
associated with the first term in Equation 5, which is 
due to rotor 1, is developed in the following. The po- 
sition vector from attachment point 1 to attachment 
point 2, expressed in the rotating hub axes, is given 
as, 


dm.i 


(*R1 “ *«) + 

rp(l) rp{ 2 ) -( 3 ) 

1 RX.i 1 Rl,i r R\ ~r 

jO) 7^(2) 7^(3) r ( r 7 A\ 

1 Rl,i 1 Rl ) i 1 Rl,i t Rl (74) 


Sri is the position vector from the center of the hub 
to attachment point 1, expressed in hub axes. t Rl 


f R\ , i = FRi ti (vRi ti ) (77) 


The hub axis components of the forces acting at the 
attachment points are given as, 


j?D 

— -(Fri ,»•/ | \)d R1>i 

(78) 

p d 

2 

= (Cri,./ I d RX% i |)rfju |t - 

(79) 

The velocities at the attachment points, expressed 

in the hub axes, 

are: 


= 

Psi7ffi,«] T (vs + 5(wf’)rf 1 ) + 



Pffi,i] T S(wsi)r/ri + 



S(uH\,i)s R \ 

(80) 

v D - 

V Rl,i, 2 — 

V Rl,i,l + 

(81) 


The above two equations are used to calculate the 
required partial velocities, and . 

Generalized Spring-Damper Forces Due to 
Torsional Spring-Dampers 
The torsional springs and dampers mounted on the 
blade hinges are assumed to possess linear stiffness 
and damping properties. They give rise to the third 
and fourth terms in Equation 5. The third term in 
this equation is due to rotor 1 and is discussed below. 
Referring to Figure 2, the following quantities are de- 
fined for the ith blade. Mr\^i denotes the body axis 
components of the moment on acting link 1 due to the 
spring-damper at hinge 1. Mr \^2 denotes the body 
axis components of the moment on acting link 2 due 
to the spring-damper at hinge 2. Mr 1|i(3 denotes the 
body axis components of the moment on acting on 
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the blade due to the spring-damper at hinge 3. These 
moments can be expressed as, 


MriS.i 

Mfll.i, 3 


= -b 
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,) (83) 

.i)(84) 


where fcp and ko denote sttifness and damping con- 
stants. The torsional spring-damper moments acting 
on the hub, link 1, link 2, and the blade, expressed in 
their individual body axes, are respectively given as, 


Mg Utl 

= — AfRl.t.1 


(85) 

MRu, 2 

= Afrti,«\i - 

2 

( 86 ) 

3 

= ^Rl,i,2 ~ 

M HI ,i,3 

(87) 


= 3 


( 88 ) 


The individual body axis components of the angu- 
lar velocities of the hub, link 1, link 2, and the blade, 
are respectively given as, 


, P 

w /ll,i,l 

= WHl.i 

(89) 

, ,D 

W H1,»,2 

- 

— w /ll,i 

(90) 

, P 

- t J 2 ) 

~ u Rl,i 

(91) 

, P 

w /ll ,i,4 

= Wfll.i 

(92) 


The preceding four equations are used to compute 
the four partial velocities needed for evaluating the 
third term in Equation 5. 


Airframe Kinematics 


To complete the description of the airframe dy- 
namic model, the kinematic relationship between the 
vectors q , q, and u needs to be stipulated. Let the 
airframe kinematic equations be given as, 


/*,(?>?>«) = 0, i = l,...,n RB (93) 


The elements of the vector fx i are given in detail as 
follows: 


\ fK : 1 
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(95) 


The matrix W B is given as, 


W B 


10 — sin 6 

0 cos <f> sin 4> cos 9 

0 — sin cos <j> cos 9 


( 96 ) 


Multiblade Coordinate Transformation 

The airframe dynamic and kinematic models, given 
by Equations 1 and 93, respectively, are derived in 
the rotating system, with the rotor degrees of free- 
dom describing the motion of individual rotor blades. 
However, the rotor usually responds as a whole to ex- 
citation and for physical insight it is desirable to work 
with the degrees of freedom which model the entire 
rotor system rather than the individual blades. To 
transform the equations of motion with respect to 
individual blade coordinates to rotor system coordi- 
nates, the method of multiblade coordinates is used 
[11]. Considering the example of rotor 1, for like de- 
grees of freedom, the kth (k = 1, . . npi) individual 
rotor blade degree of freedom is expressed as, 

<■ = am,o+ 


(n«i-l)/2 

y; (“fli.ic cos iipk + ciRi'i, sin iif> k ) (97) 

i=i 

for rotors with an odd number of blades and 
= ap^o-b 


(n/ji-2)/2 

c cos ixfok +aR\ t is sin iifik) + 

1=1 

(98) 

for rotors with an even number of rotor blades. 

Let the generalized coordinate and generalized 
speed vectors in the multiblade or non-rotating co- 
ordinate system be represented by q' and u f . Then 
the following substitutions are made in the airframe 
kinematic and dynamic model descriptions, given by 
Equations 93 and 1, respectively. 


T(*W 

(99) 

T(t)q' + T(t)u' 

(100) 

T{t)q' + 2T(t)u' + T{t)u' 

(101) 


Then the resulting airframe kinematic and dynamic 
equations can be written as, 

w') = 0, i = 1, . . .,n RB (102) 
/.(g'.u'.iu.u',*) = 0, i = 1 ,... ,n RB (103) 
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where the vector w consists of the inflow coordinates 
of the two rotors: 

w = WRl,ls, WRl,lc ^i?2,0> ^R2,1 c} T 

(104) 


Rotor Dynamic Inflow Model 

The rotor dynamic inflow model used in this work 
is based on the Peters and HaQuang [12] model which 
is in turn based on the work of Pitt and Peters [8]. 
The model includes three inflow degrees of freedom 
that yield the time-varying induced flow parallel to 
the rotor shaft. Based on the small perturbation 
potential flow equations, the model accounts for dy- 
namic changes in collective inflow and first harmonic 
inflow azimuthally. Inflow along the blades varies lin- 
early. The inflow distribution is given by Equation 
50. For simplicity only the dynamic inflow model for 
rotor 1 will be described. The dynamic inflow model 
for rotor 2 is similar with obvious changes. 

The basic model formulation is carried out in the 
rotor wind axis system and is later transformed to the 
rotor shaft axis system. The dynamic inflow equa- 
tions are forced by the averaged (over rotor revolu- 
tion) rotor thrust, rolling moment, and pitching mo- 
ment in the shaft axes. The resulting equations can 
be written in the form, 


f WRI,0 


™Rl,0 


-Tr, ’j 
7^r f 

| WRlAs 

[ V>Rl t lc J 
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pint | 

M/ti 

l J 


(105) 


where, 

Ar\ = ARi(q',u f ) (106) 

The blade element forces, given by Equation 60, are 
vectorially summed over all rotor blades to obtain 
the shaft axis components of the rotor thrust, rolling 
moment and pitching moment. These forces and mo- 
ments are then averaged over the period of revolution 
of the rotor and used in Equation 105. 

The complete set of dynamic inflow equations for 
the two rotors can be functionally represented as, 

9i(q',u’,w,w) = 0, i= (107) 

where nm — 6 since three state inflow models are 
being used for each rotor. It should be noted that 
Equation 107 is written using the multiblade or non- 
rotating coordinate system. 


Trim Algorithm 


Trim of an aircraft is defined as an equilibrium 
condition where the translational and rotational ac- 
celerations of the fuselage are zero. Hence in trim, 
p = q= r = u = v = w = 0. For straight and 
level flight, p=g = r = v = 0as well. For a fixed 
wing airplane this definition is sufficient since one can 
generally regard an airplane as a single rigid body 
with six degrees of freedom. For rotorcraft the con- 
cept of trim is more complicated because the vehicle 
is represented as a multibody system consisting of a 
fuselage, many rotor blades, and a drive system. By 
virtue of the rotor rotational motion, the blades are 
always accelerating. For the rotor blade degrees of 
freedom, trim is considered to be an operating condi- 
tion such that the individual rotor blades follow a pe- 
riodic path. This implies that all the first and second 
derivatives of the rotor multiblade coordinates must 
be zero in trim. This will force individual blades to 
track the same periodic path each rotor revolution. 
However, it should be noted that for even bladed ro- 
tors this condition will not force every blade on a 
rotor to follow the same path. This is due to the 
warping multiblade coordinate mode for even bladed 
rotors. 

There are many different methods for obtaining 
the trim condition of a coupled rotor and fuselage 
combination. Included in these methods are itera- 
tive fuselage trim and rotor trim, fully coupled au- 
topilot trim, finite elements in time trim, nonlinear 
optimization trim, and Galerkin method trim. While 
no one method for trim is superior in all settings, all 
the methods are sufficiently different to have qualities 
which make them more or less attractive in different 
settings. In this work a nonlinear optimization trim 
technique is used. 

In nonlinear optimization, one seeks to minimize or 
maximize a certain nonlinear function by iterating on 
the independent variables of the problem. Here the 
sum of the squares of the dynamic equation residu- 
als will be minimized and the independent variables 
will be the system states and controls. A modified 
Newton’s method, sometimes called a damped New- 
ton’s method or a quasi Newton method, is used as 
the nonlinear optimization algorithm to compute the 
trim state of the vehicle. 

The trim algorithm begins by noting that in trim, 
u* = 0 and w = 0 necessarily. Hence in trim, the 
airframe dynamic equations (Equation 1) and the dy- 
namic inflow equations (Equation 107) can be written 
as, 


/»0M) 
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0, i = 1, . . . , n RB (108) 

0, * “ !,*•■» n Di (109) 
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where 


x = {(< ? ') T ,(u / ) T .H T } T (110) 

Clearly, x is the state vector of the airframe dynamic 
model. Equation 108 contains a set of algebraic non- 
linear equations which are periodic in time, with a 
period of r, r is the period of revolution common to 
rotor 1 blades and rotor 2 blades. The goal of the 
trim algorithm is to minimize the residual of each 
equation in Equations 108 and 109 for all values of 
time. 

A natural scalar function to minimize for trim is, 


1 r T Ti RR up/ 

J = - I ^ + 

T J ° . = I >=t 
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/L, -C >(^) 2 + (Hi) 


Jfc = l »=sl 


*=1 


where nj is the number of time points chosen for dis- 
cretization. The function J is termed the cost func- 
tion. Using the discretized form of the cost function, 
the gradient and hessian of the cost function can be 
formed. 
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The minimization problem described above is essen- 
tially a least squares problem. It is known that for 
least square minimization problems, where the cost 
function is small at the solution, the second deriva- 
tive terms in the above equations are relatively small 
and can be neglected [13]. By definition, this assump- 
tion is valid in the trim problem. 

In a modified Newton’s method, a local optimiza- 
tion problem is solved iteratively. A flow chart for 
the iteration procedure is given in Figure 6. Using 
an initial condition or guess for the trim variables, a 
local quadratic model of the cost function is formed, 


J(x + Ax) = J(x) + ^Ax + (114) 

At the local minimum of this approximation to the 
actual cost function one must have, 


dJ 

dAx 


= 0 


(115) 


For a local minimum of a quadratic function to exist, 
hessian matrix of the cost function must be positive 
definite. Assuming this is the case, 


Ax = — 


' d 2 J ' 
dx 2 



(116) 


The vector Ax is called the search direction because 
based on this direction a search to reduce the cost 
function shall be undertaken. For the local quadratic 
model of the cost function, the minimum is given by 
x+ Ax, of course if a minimum exists. A new iteration 
on the minimum of the actual cost function can now 
be made by with the equation, 


£new — %old *"f“ Oc Ax (117) 

The parameter, a, is the step length. It is used be- 
cause the local model is only an approximation to 
the actual cost, or = 1 corresponds to a full Newton’s 
method while a < 1 implies a damped or modified 
Newton’s method. The parameter a is determined at 
each trim iteration and is based on satisfying criteria 
for tracking sufficient decrease in the cost function at 
each iteration in the overall minimization problem. 
The process of determining the step length is called 
a step length procedure or line search strategy. 

There are many criteria for determining suffi- 
cient decrease in the cost function at each iteration. 
Armijo’s rule is used here which can be stated as, 

Q j 

Jo — Ja > — /Ja— Ax (118) 

where the constant /i is a positive number. A back 
tracking strategy is used in the line search strategy. 
In this method, one always starts with a — 1 and 
tries to use the full Newton’s method if possible. If 
the current a does not fulfill the Armijo condition, 
then a is divided by a factor and retried. Once an 
appropriate value for a is obtained, new values for 
x are computed. Then a new local quadratic mode! 
is formed and the optimization procedure is again 
formed. It should be noted that in solving for the 
search direction a linear system must be solved. It is 
solved using a modified Choleski decomposition algo- 
rithm as described in reference [13]. 


Linear Model of Airframe Dynamics 

Linearized rotorcraft dynamic models are ex- 
tremely useful for flying qualities analyses. To this 
end, the composite airframe dynamic model consist- 
ing of the kinematic, dynamic, and dynamic inflow 
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models, given by Equations 102, 103, and 107, re- 
spectively, is linearized about an arbitrary trim state, 
xq. The linear model can be written as, 


C p (xoJ)6x = D p (xo ) t)6x (119) 


The (2 tire + noi) x (2 urq -f rtoi) square matrices 
C p and D p are given as, 


C P 


Dp 


where 


Sk 

f 

9 


W o o 


0 

0 

1 


(120) 

r mjl 

dq' 

A 

dq* 

*JL 

L dq 1 

1 1 

(121) 

{fKx , ■ 

• > fK nRB } 

r 

(122) 

{/I.--- 

f \ T 

> JnRB J 


(123) 

{?!,••• 

\T 

j 9noi J 


(124) 


In the ensuing analysis, the <5’s in Equation 119 will 
be dropped and the perturbation state of the aircraft 
will be simply denoted as x ac . 


Transformation of the Airframe Linear 
Dynamic Equations 

The multiblade coordinate transformation should 
be accompanied by a transformation of the equa- 
tions of motion to the non-rotating coordinate sys- 
tem. This step is accomplished by taking linear com- 
binations of the equations of motion given by Equa- 
tion 119. The operations can be performed by pre- 
multiplying the dynamic equations by a transforma- 
tion matrix, T(t). The fully transformed linear equa- 
tions are, 


T(t)C p (t)x ac = T(t)D p (t)x ac (125) 

In rotorcraft handling qualities analysis, a linear 
time invariant system is most convenient to work with 
due to the powerful linear system analysis tools avail- 
able. A standard approximation used in rotorcraft 
handling qualities work is to neglect the harmonic 
content in Equation 125 and hence obtain a linear 
time invariant system. This approximation is known 
as the constant coefficient approximation and it is 
used in the current effort. 

The blade pitch control terms can be separated 
from the above equations by assuming that the multi- 
blade coordinate blade pitch degrees of freedom do 


not possess dynamics. Appropriate rows of the dy- 
namics matrix are deleted and the associated columns 
form the controls matrix. The final form of the air- 
frame linear dynamic equations is, 

x ac = Ax ac + Bd (126) 

Vac = CXac+Dti (127) 

where the vector d consists of individual rotor pitch 
control variables. This system can now be coupled to 
the flight control system to form the complete system. 


Linear Control System Model 

Most aircraft flight control systems are given in 
block diagram form and there is no standard struc- 
ture. Although for modeling purposes, a generic flight 
control system structure could be assumed such that 
all or at least a majority of current aircraft flight con- 
trol systems could be accommodated, it is felt this 
approach may be too restrictive in some cases and 
far too general, hence inefficient, in other cases. It is 
desirable to have a flight control system modeling ca- 
pability which does not assume a structure aprior but 
uses the input data deck to generate the model. This 
approach allows for greater flexibility and increased 
utility of the control system model. With these con- 
siderations in mind, a linear state-space flight control 
system modeling capability was developed that takes 
the basic block diagram data as input. 

The flight control system is assumed to be com- 
prised of an arbitrary number of filters, given in poly- 
nomial form. Each filter is a multi input and single 
output filter as shown in Figure 7. 

The inputs to each filter can consist of pilot stick 
inputs, outputs of other individual filters, aircraft 
states, and derivatives of aircraft states. A state- 
space realization is computed for each individual fil- 
ter in phase variable canonical form. The filters are 
then assembled into an overll state-space realization. 
The realization can be written as, 

ics + B u Vu (128) 

y u = C u x cs + D u v u (129) 

The subscript u signifies that the state-space matri- 
ces do not account for the filter coupling. A filter 
coupling matrix can be computed in the form, 

— CuVu "b ~b ^u^ac (130) 

It should be noted that Equations 128, 129 and 130 
can be constructed in a straight forward manner from 
the input block diagram data. Substituting Equation 
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130 into Equations 128 and 129, the coupled state- 
space model of the control system can be formed. 
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The matrix X restricts the overall control system out- 
puts to be the aircraft blade pitch angles. It should 
be noted that if the matrix [/— £> u ( u ] is singular, then 
there is not a valid state-space model for the system 
and the system is non-causal. This is due to the fact 
that the flight control system output can be written 
as, 

A.C«] V — G% i%cs 4 &uftu & 4 Du 4 Du&u Xqc 

(I 45 ) 

For a valid state-space realization the output must 
be uniquely determined from the state and control. 
Clearly when [/— D u £ u ] is singular this is not possible. 
This observation can be used for detecting input data 
errors. 

Airframe and Control System Coupling 

The linear airframe model which describes the rigid 
body aircraft motion and rotor dynamic inflow is 
given by Equations 126 and 127. The inputs to the 
airframe linear equations are blade pitch angles of 
the two rotor systems. The linear flight control sys- 
tem model is given by Equations 131 and 132. The 
outputs of the flight control system model are also 
the blade pitch angles of the two rotors. The linear 
airframe and control system models are coupled by 
noting that the output of the flight control system 
model is the input to the airframe model. 
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where, 



An = 

a + b n 

(148) 

A 12 = 

ST 

(149) 

A 21 = 

H + E(A + B n) 
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(151) 
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Ci = 

c + d n 

(154) 

c 2 = 

D T 

(155) 

Dx = 

DE 

(156) 

n = 

| [/■ -ZB]' ~ l (R- +ZA) 

(157) 

T = 

[I - ZB]~ l P 

(158) 

E = 

[I — ZB]~ l Q 

(159) 


Concluding Remarks 

A linear coupled rotor-fuselage-control system dy- 
namic model is presented in this paper. The model is 
expected to be useful for flying qualities studies, sta- 
bility and control investigations, and control design 
parametric studies. Efforts are underway to produce 
numerical results for the validation of the model. 
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The matrices E x , £ 2 , and £ 3 represent single axis 
transformations about x, y , and z axes, respectively, 
and are defined as follows: 


E,{k) 


£ 2 («) 


E 3 (k) 


10 o' 

0 cos k sin k 
0 — sin k cos k 

cos k 0 — sin k 

0 1 0 

sin k 0 cos k 

cos k sin k 0 

— sin k cos k 0 

0 0 1 


Time-Derivative of a Transformation Matrix 
Consider the time-derivative of a vector v in two 
reference frames denoted by A and B. Let the com- 
ponents of v in Frame A be denoted by v A and those 
in Frame B be denoted by vq, Let the angular ve- 
locity of Frame B with respect to Frame 4 be w and 
let the components of u in Frame B be denoted by 
Wfl. Let T represent the transformation matrix that 
transforms vector components from Frame B axes to 
components in Frame A axes. The time-derivatives 
of v in Frame A and Frame B are related by the fol- 
lowing vectorial equation: 


A dv 

dt 


R dv 

= Tt +U * v 


In matrix-vector format, the preceding equation can 
be written as, 


v A _ Tvq - f - TS{u>b)vb 

Also, since v A = Tvq , one gets for i ) A the following 
expression: 

v A = Tvb + Tvq 

Comparing the two equations for v A) the following 
formula is obtained for T : 


T = TS(u> b ) 


Appendix 

Skew-Symmetric Matrix Operator 

For a vector a = {ai , a 2 , a 3 } T , the matrix 5(a) is 
defined as, 


S(a) 


0 -a 3 a 2 

a 3 0 -ai 

— a 2 a x 0 


Single Axis Transformation Matrices 
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Figure 2: Rotor Blade Geometry 


2 



Figure 6: Trim Procedure Flowchart 



y (1) 


ni(i) « Number of Inputs to i'th Filter 
no (i) - Order of i'th Filter 
yU) - Output of i’th Filter 

- Gain Multiplying j'th input of i’th Filter 

Figure 7: zth Flight Control System Filter 
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